#This to obtain:
#table_A12.tex


library("broom")
library("dplyr")
library("grid")
library("tidyr")
library("sandwich")
library("lfe")
library("broom")
library("dplyr")
library("tidyr")
library("grid")
library("stargazer")
library("lmtest")
library("multiwayvcov")
library("ggplot2")
library("glue")
library("spdep")
library("spatialreg")
library('dismo')
library("stargazer")
library("rgdal")
library("readxl")
library("stringr")
library(spatialEco)

#setwd("C:/Users/bogdanp/")
setwd("/Users/bgpopescu/")


#Step 1: Read Data
census_1950_dist <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                     layer="croatia_district_centroids_1950_GCS")
plot(census_1950_dist)


#Step2: Making polygons
krajna6_treat_HABS <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                              layer="krajna6_treat_HABS_GCS")
krajna6_treat_HABS<-subset(krajna6_treat_HABS, select = c(NAME_0))

#Step3: Making function for geometry
rename_geometry <- function(g, name){
  current = attr(g, "sf_column")
  names(g)[names(g)==current] = name
  st_geometry(g)=name
  g
}

#Step4: Making another polygon
krajna6_treat_HABS = rename_geometry(krajna6_treat_HABS, "SHAPE")

krajna6_treat_KRAJ <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                              layer="krajna_polygon_d_GCS")
krajna6_treat_KRAJ<-subset(krajna6_treat_KRAJ, select = c(name_regiment))
st_geometry(krajna6_treat_HABS) <- "SHAPE"
names(krajna6_treat_KRAJ)[names(krajna6_treat_KRAJ)=="name_regiment"]<-"NAME_0"
krajna6_treat_KRAJ = rename_geometry(krajna6_treat_KRAJ, "SHAPE")

#Step5: Merging the polygons
HRV_adm_comp <- st_combine(rbind(krajna6_treat_HABS, krajna6_treat_KRAJ))
HRV_adm_comp<-st_make_valid(HRV_adm_comp)

#Step6: Creating Voronoi Polygons function
st_voronoi_point <- function(points){
  ## points must be POINT geometry
  if(!all(st_geometry_type(points) == "POINT")){
    stop("Input not  POINT geometries")
  }
  g = st_combine(st_geometry(points)) # make multipoint
  v = st_voronoi(g)
  v = st_collection_extract(v)
  return(v[unlist(st_intersects(p, v))])
}


#Step7: Using the Voronoi Polygons function
p<-st_as_sf(census_1950_dist, remove = FALSE)
v = st_voronoi_point(p)
pv = st_set_geometry(p, v)
x <- st_intersection(st_cast(pv), st_union(HRV_adm_comp))
plot(x)

census_1950<-read_excel("./Dropbox/Legacies_Project/Analysis/data/census_1950.xlsx",
                     sheet=1, col_names = TRUE, skip = 0)

census_1950 <- subset(census_1950, country == "Croatia")
census_1950<- merge(x, census_1950, by.x="Name_1950", by.y="name")


census_1950$treat<-census_1950$inside_krajna6_treat
census_1950$dist1 <- as.numeric( with (census_1950,ifelse(census_1950$treat==1, 1, -1)))
census_1950$krajna6_distance<-census_1950$krajna6_distance/1000
census_1950$zagreb_distance<-census_1950$zagreb_distance/1000

census_1950 <- subset(census_1950, inside_krajna6_ctrl == "1" | inside_krajna6_treat == "1")


census_1950$dist2 <-census_1950$dist1*census_1950$krajna6_distance
census_1950$bfe1 <- ifelse(census_1950$krajna6_NEAR_FID == 1, 1,0)
census_1950$bfe2 <- ifelse(census_1950$krajna6_NEAR_FID == 3, 1,0)
census_1950$bfe3 <- ifelse(census_1950$krajna6_NEAR_FID == 4, 1,0)
census_1950$bfe4 <- ifelse(census_1950$krajna6_NEAR_FID == 5, 1,0)
census_1950$bfe5 <- ifelse(census_1950$krajna6_NEAR_FID == 6, 1,0)
census_1950$bfe6 <- ifelse(census_1950$krajna6_NEAR_FID == 7, 1,0)
census_1950$bfe7 <- ifelse(census_1950$krajna6_NEAR_FID == 8, 1,0)
census_1950$bfe8 <- ifelse(census_1950$krajna6_NEAR_FID == 9, 1,0)


census_1950$OttZdist=census_1950$treat*(census_1950$zagreb_distance)
census_1950$OttZdist2=census_1950$treat*(census_1950$zagreb_distance^2)

treat_group <- census_1950[ which(census_1950$treat==1),]
control_group <- census_1950[ which(census_1950$treat==0),]



census_1950$quad1<-ifelse((census_1950$POINT_X > 18 & census_1950$POINT_X< 20) & 
                       (census_1950$POINT_Y<47  & census_1950$POINT_Y>44), 1, 0)

census_1950$quad2<-ifelse((census_1950$POINT_X > 16 & census_1950$POINT_X< 18) & 
                       (census_1950$POINT_Y<47  & census_1950$POINT_Y>44), 1, 0)

census_1950$quad3<-ifelse((census_1950$POINT_X > 14 & census_1950$POINT_X< 16) & 
                       (census_1950$POINT_Y<47  & census_1950$POINT_Y>44), 1, 0)


column_labels <-c("\\shortstack{Total Cooperatives\\\\ 1950}",
                  "\\shortstack{Households with No\\\\ Cooperative Members\\\\ 1950}",
                  "\\shortstack{Cooperative Plows\\\\ 1950}")


zadrugas_total<-lm(zadrugas_total~treat + POINT_X + POINT_Y + zagreb_distance + 
  bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7+
  quad1 + quad2 + quad3, data=census_1950)
summary(zadrugas_total)

households_with_no_zadruga_member<-lm(households_with_no_zadruga_member~treat + POINT_X + POINT_Y + zagreb_distance + 
                     bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7+
                     quad1 + quad2 + quad3, data=census_1950)
summary(households_with_no_zadruga_member)

coop_machines_plows_all_kinds<-lm(coop_machines_plows_all_kinds~treat + POINT_X + POINT_Y + zagreb_distance + 
                                        bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7+
                                        quad1 + quad2 + quad3, data=census_1950)
summary(coop_machines_plows_all_kinds)

communist =stargazer(zadrugas_total,
                     households_with_no_zadruga_member,
                     coop_machines_plows_all_kinds,
                     column.labels= column_labels,
                     keep = c("treat"),
                     covariate.labels=c("Military area"),
                     se= list(coeftest(zadrugas_total, cluster.vcov(zadrugas_total, census_1950$Name_1950))[, 2],
                              coeftest(households_with_no_zadruga_member, cluster.vcov(households_with_no_zadruga_member, census_1950$Name_1950))[, 2],
                              coeftest(coop_machines_plows_all_kinds, cluster.vcov(coop_machines_plows_all_kinds, census_1950$Name_1950))[, 2]),
                      dep.var.caption = "Dependent variable:",
                      dep.var.labels.include = F,
                      omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                      add.lines=list(c("Mean", round(mean(census_1950$zadrugas_total, na.rm=T), digits = 3),
                                       round(mean(census_1950$households_with_no_zadruga_member, na.rm=T), digits = 3),
                                       round(mean(census_1950$coop_machines_plows_all_kinds, na.rm=T), digits = 3)),
                                     c("SD", round(sd(census_1950$zadrugas_total, na.rm=T), digits = 3),
                                       round(sd(census_1950$households_with_no_zadruga_member, na.rm=T), digits = 3),
                                       round(sd(census_1950$coop_machines_plows_all_kinds, na.rm=T), digits = 3)),
                                     c("Boundary FE", rep("Yes", 7)),
                                     c("Region FE", rep("Yes", 7))))
note_latex <- "\\multicolumn{4}{l} {\\parbox[t]{16cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and robust standard errors in parantheses from OLS regressions. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The unit of analysis is canton (srez). 
See codebook for data sources and how the variables were calculated.}}"
communist [grepl("Note", communist)] <- note_latex
cat(communist, file="./Dropbox/Legacies_Project/Paper/tables/table_A12.tex", sep="\n")
